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DISCUSSION PAPER 
CONDITIONAL GROWTH CHARTS^ 

By Ying Wei and Xuming He 

Columbia University and University of Illinois 

Growth charts are often more informative when they are cus- 
tomized per subject, taking into account prior measurements and 
possibly other covariates of the subject. We study a global semipara- 
metric quantile regression model that has the ability to estimate con- 
ditional quantiles without the usual distributional assumptions. The 
model can be estimated from longitudinal reference data with irreg- 
ular measurement times and with some level of robustness against 
outliers, and it is also flexible for including covariate information. We 
propose a rank score test for large sample inference on covariates, and 
develop a new model assessment tool for longitudinal growth data. 
Our research indicates that the global model has the potential to be 
a very useful tool in conditional growth chart analysis. 

1. Introduction. Growth charts, also known as reference centile charts, 
have been widely used to screen the measurements from an individual sub- 
ject in the context of population values. A growth chart consists of a series 
of smooth curves plotted against time or another covariate, with each curve 
representing the trend of a given percentile of the measurement in a pop- 
ulation. With several chosen percentile curves, a growth chart displays the 
distribution of a certain measurement within a certain range of time for a 
certain population. When a measurement is extreme on the chart, the sub- 
ject is often identified for further investigation. An extreme measurement is 
likely to be a reflection of some unusual underlying physical condition. 
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The conventional method of constructing reference centiles is to get the 
empirical percentiles from cross-sectional data at a series of time points, 
and then fit a smooth polynomial curve to them. This method was used to 
develop the National Center of Health Statistics (NCHS) Growth Chart in 
1977; see [15]. Cole [5] has become a classical work in the statistics liter- 
ature on growth chart construction. In recent years, a number of different 
methods have been developed in the medical statistics literature; see, for 
example, [1, 9, 22, 32, 34, 39, 42]. Wright and Royston [43] have reviewed 
some of these methods. Most authors explore the beauty of normality to 
reduce quantile estimation to estimation of moments. Since most physical 
measurements are known to be nonnormal in distribution, a transformation 
to normality is generally used. The Box-Cox power transformation remains 
the most popular choice in this regard, as it can be found in a wide variety 
of medical journals — [7, 8, 29, 30], just to name a few. Arguably the most 
sophisticated transformation method is the LMS method of Cole [5]. It as- 
sumes that, at any time t, the measurement Y{t) can be transformed to be 
approximately standard normal by Z{t) = ((y(t)/M(t))'^(*) - l)/{L{t)S{t)), 
where L(t) is the best-fitting Box-Cox power for Y(t), M{t) is the median 
of Y[t)^ and S{t) is the scale of the transformed variable. The transformed 
variable Z{t) is also known as the z-score. The penalized log-likelihood ap- 
proach for estimating the L-M-S functions has been proposed by Cole and 
Green [9]. 

Since conventional growth charts are usually developed from cross-sectional 
data, they are most useful for examining a subject with one measurement 
at a specific time. If a subject has more than one measurement, it can be 
more informative to study his growth path rather than a single measure- 
ment. Knowing a subject's prior growth path gives us a better understand- 
ing about his current growth status. For example, if a subject who has been 
growing along the 75th percentile curve suddenly drops below the median 
curve, it might be worth singling him out for further investigation. If we sim- 
ply compare the growth path against the reference centiles generated from 
cross-sectional data, it would be hard to know how big a drop should be 
deemed unusual. Another important reason arguing against this approach is 
given in [6]. During both infancy and puberty, the subjects on the upper (or 
lower) centiles tend to grow toward the median at a faster growth rate than 
others, which is known as "catch-up" growth. Their growth paths are thus 
likely to move across the reference centiles during infancy or puberty. This is 
a normal phenomenon, but tracking those subjects on a conventional growth 
chart may give us an incorrect impression of abnormal growth. Clearly, one 
should compare a subject with the group of subjects with similar growth 
paths. In this paper, we focus on developing conditional growth charts based 
on longitudinal data. 
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Several authors, including Royston [33], have considered models for con- 
ditional reference centiles. Cole [6] considered the LMS-AR model Zt = 
at + btZt-i + et, where Zt is the so-called z-score from a power transfor- 
mation at time t, Zt-\ is the lagged measurement at time t — \ and ej is 
distributed as A^(0,(Tf). In addition to the assumption of normality, this 
method requires fixed measurement time intervals. A more general model 
allowing varying measurement time intervals has been proposed by Thomp- 
son and Fatti [38], but they assumed a multivariate normal distribution for 
the measurements and the covariates at all time points and used the maxi- 
mum likelihood estimator for the mean and variance functions. Scheike and 
Zhang [35] and Scheike, Zhang and Juul [36] considered a longitudinal re- 
gression model = m(lij_i, tjj — iij-i) for height or weight, where 
lij is the jth measurement of the ith subject at a random time tjj-, m(-) 
is an unknown function of the prior measurement and the time duration 
from the prior measurement, and e^j- is normally distributed. The reference 
centiles are constructed based on an estimate of m and the normality of 1^ j- . 

In practice, the longitudinal data collected for constructing reference cen- 
tiles often exhibit the following features. First, the measurement time in- 
tervals are somewhat random. Even under the standard guideline of taking 
measurements on a nationally accepted schedule, say bimonthly for young 
children, the actual measurement times deviate from the fixed schedule due 
to practical considerations. A new subject to be screened for medical reasons 
is unlikely to follow the same schedule as the reference group. Second, trans- 
formation to normality is often a reasonable endeavor at a given time point, 
but the conditional models would require more than marginal normality of 
the z-scores. Normality of the error distributions in the conditional models 
can fail even if the marginal measurements are nearly normal; see an example 
later in the paper. Third, it is often desirable to incorporate subject-level co- 
variates such as parent height in the conditional growth analysis. Normality 
in the conditional models becomes a stronger requirement in the presence of 
covariates. The existing approaches to conditional reference centiles can be 
ineffective in view of the data features discussed here. The semiparametric 
quantile regression model proposed in [40] aims to better accommodate data 
of this type; see illustrations in [41]. 

In this paper we further develop the global semiparametric approach to 
conditional growth charts by providing asymptotic theory, inference tools 
and a model assessment technique. Section 2 describes the global model and 
a simple estimation procedure. Section 3 provides the large sample properties 
of the estimates, and Section 4 gives a rank score test for covariate effects in 
the model. Section 5 compares the global model with the LMS-AR method of 
Cole [6] and illustrates the difficulty in preserving normality in conditional 
models. We also provide a new tool of model assessment for longitudinal 
growth data in Section 6. Some practical issues of using growth charts, both 
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conditional and unconditional, are given in Section 7, and the technical 
proofs of our theorems are given in Section 8. 

2. A global model. Suppose that we have n subjects, and the ith. subject 
has rrii measurements at times ii,i,ti,25 • • • j^i.mii which are not necessarily 
evenly spaced. We denote as Yij the jth measurement of the ith subject, 
and Di^j^k = ti,j — ii,j-k is the time distance between the jth and {j — k)th 
measurements. For a given r, we model Yij by 

p 

(2.1) Y,j = goiUj) + Y.{ak + 6fc A,i,fc)>i,i-fc + Xj^-fo + eijir), 

k=l 

where i = 1, . . . , n, j = p + 1, . . . , m^, Xij = (X^j^i, . . . , Xij^i^ consists of / 
covariates for the ith subject at time Uj, and eiji^r) is a random variable 
whose rth quantile, given the growth path up to the {j — l)st measure- 
ment and Xij, is zero. The conditional quantile of Yij given the p prior 
measurements and the covariate Xij includes go as a nonparametric inter- 
cept function of the current measurement time, an autoregressive function 
of Yij-i, . . . ,Yij-p whose coefficients are linear functions of measurement 
time distances Dij^k, and a linear function of the covariate Xij. 

In general, the subjects take measurements at irregular time intervals. 
It is natural to assume that the dependence between two measurements 
varies with their measurement time distance. This motivates the choice of 
the autoregressive coefficients as functions of Di j^j.. If all the subjects take 
measurements at the fixed time intervals, the coefficients simplify to con- 
stants. On the other hand, it is possible to generalize the model by using a 
more general form as considered by Cai and Xu [2] . We defer this discussion 
to Section 7. 

We have intentionally written the error term in model (2.1) as ejj(r), 
which is r-specific. This helps distinguish the quantile regression model from 
most other models of the same form. Model (2.1) allows the error term to be 
dependent on the covariates as well, so eij{T) is merely a convenient notation 
for the difference between Yij and its rth conditional quantile function. If 
the error terms are independent of r and of the covariates in the model, 
then models of this form have been well studied in the longitudinal data 
analysis by many other authors, including Chiang, Rice and Wu [4], Lin and 
Carroll [24], Lin and Ying [23] and Fan and Li [11], to mention just a few. 

For the sake of simpler presentation in the rest of the paper, we will 
convert to the usual notation of writing eij(r) as eij, but it is always helpful 
to keep in mind that this error term is r-specific. 

In (2.1), go{t) is a nonparametric function with a certain degree of 
smoothness. Suppose that the range of interest in time is t € [tL,tu]- Dif- 
ferent smoothing methods may be used to estimate go, but Wei [40] used 
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the convenient approach of regression splines. Specifically, we approximate 
go{t) by a linear combination of B-spline basis functions as in [18]. Given 
a set of knots, we denote as 7r(t) = (7ri(t),7r2(t), . . . ,7rfc^(t))''' the set of kn 
B-spline basis functions of, say, order 4 (cubic). Let TTjj- = 7r{tij), Xij = 
{Yij-i,Dij^iYij-i, . . . ,Yij-p,Dij^pYij-p,Xi^jy; we can then use nj^an + 
XjjjS^ to estimate the rth conditional quantile function of Yi j, where q;„ 
and /3„ are obtained by minimizing 



where Prir) = rir — I{r < 0)) is the quantile objective function; we refer 
to [26] for details. 

Obviously, other smoothing methods (e.g., kernel smoothing) may be used 
in lieu of regression splines. One advantage of spline smoothing here is that 
in typical applications it is usually sufficient to preselect a set of knots in 
the interval (tL,tu) using our general understanding of growth patterns. 
For example, it would be useful to place more knots during infancy and 
puberty than during other times. In this paper, we do not go into the issue 
of automated knot selection. 

Optimization of (2.2) can be performed efficiently by linear programming 
techniques. In fact, the solution to (2.2) only requires software for quantile 
regression in linear models. The R package quantreg contributed by Koenker 
and the SAS PROG QUANTREG are handy for obtaining the estimates of 
a and (3. 

We further note that in model (2.1), even if the e^j's are assumed to be 
independent, the within-subject correlation may be captured by the autore- 
gressive component of the model. For a consistent estimate of the conditional 
quantile function, we do not need the i.i.d. assumption on e^j. The inference 
tools developed later in this paper are proven under stronger assumptions 
but with good robustness against common deviations from the i.i.d. assump- 
tion. 

3. Large sample property. Under appropriate conditions, the intercept 
function estimate and the coefficient estimate for (3 converge to their true 
values at the optimal rates. To state the conditions, we use || • || for the 
Euclidean norm, and ctn ~ to mean < liminf^ af),/6n, < limsup^j (Xn,/ftn ^ 
cxD. Throughout, we assume that the number of measurements per subject 
is bounded, but the number of subjects grows. For the sake of convenience 
and without loss of generality, we assume that p = 1 in model (2.1), and 
consequently drop the subscript k from Dij^k so that -Djj = Dij^i. Also let 
'>Pix) = p't{x) = T-/{i.<o}, and -(/-(e,) = (^^(64,2), ^'(ej.a), • • • , ^(ei.mj)'^- Since 



(2.2) 




i=l j=p+l 
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our results will be stated for a given r, we have dropped the dependence of 
tp on T here. 

Following He and Shi [18], we assume a spline approximation to 50 (i) as 
TTj^jQoi with Rnij = T^JjOLQ — gQ{tij) as the approximation error. We state 
the assumptions for this section. 

Assumption 1. goii) has bounded rth derivative for some r > 1. 

Assumption 2. The numbers of measurements rrii are uniformly bounded 
for all n. 

Assumption 3. The errors Cj = {ei^2, ■ ■ ■ , Gi,mi)'^ have unspecified variance- 
covariance structures but are independent across subjects. 

Assumption 4. Ai = E[tp{ei)ip{eiy] > for each i and supj ||yli|| < 00. 

Let bij = /eij (O), where /e^^ is the conditional p.d.f. of eij given tj.j and 

Assumption 5. < infijbij < sup^ jbij < 00, and as s ^ 0, 
supi j \E'ip{eij + s) - bijs\ = O(s^). 

Assumption 6. There exists an (/ + 2)-dimensional function r]{t) = 
(771 (i), 1^2(0) ■■■) ^/+2(i)) with bounded rth derivative such that for any i 
and j, 

where 6ij = (5ij,2, ■ • ■ , Si,j,i+2) are random vectors with mean zero and 

are independent of e^j given Let Xn = (-^1,2, • ■ • ,-^i,mi,-^2,2, • ■ • ,-^n,m„) 
be the {I + 2) x N matrix with N = J2i=i{'mi - 1)) and let A„ and K be 
defined in a similar way using 6ij and 'r]{tij), respectively. Then we have 
= Ki + Am where is the mean of so that EAn = 0. 

Assumption 7. sup„n~^£'||A„|p < cx). 

Assumption 7 is readily satisfied if all the variables have finite supports. 
Note that for any i, (ti,i,ti,2; ■ • • ^U,mi) are order statistics of the time vari- 
able, thus tjj-i and Yi,]-! are not independent of tjj. Therefore, rji{ti^j) 
is the conditional mean of the prior measurement given the current time, 
and T]2{tij) is the conditional mean of Dij times Yij-i given current time 
tjj. Similarly, (?73, . . . ,7?/+2) are the conditional means of the covariates 
given tij. 
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Let Zn — (vri,2,7ri^3, . . . ,7Tij, . . . ,7rn,m„)7VxA;„' ^ri — ^n(^n ^n)^n ^^^'^ — 

diag(6ij), where {hj} is sequenced by indices i and j as in Zn- For additional 
assumptions, we define X* = {I — Gn)Xn and Kn = X*^ BnX*. If we denote 
as X*^'^ the nii X (I + 2) submatrix of X* corresponding to the zth subject, 

and = diag(6,,2, • • • , then K„ = E?=i Xn^'^^S^^)^*^') . By the as- 

sumption of between-subject independence, Kn can be regarded as an inde- 
pendent sum. Furthermore, let A„ = diag(^i, . . . ,An) and Sn = X*^ AnX*. 
Then we have 

Assumption 8. There exists a positive definite matrix K such that 
Kn/n K in probability, and E{Kn/n) K. 

Assumption 9. There exists a nonnegative definite matrix S such that 
Sn/n^ S in probability. 

Let gn{U,j) =T^J^jO'-n\ the following theorem summarizes the large sample 
properties of /3„ and 

Theorem 3.1. (i) Under Assumptions 1-7, if the number of knots kn ~ 
^i/(2r+i)^ i/ien 

(3.1) ||^„-/3o||=Op(n-'-/(2'-+i)) 
and 

1 n nii 

(3-2) n^T. i9niU,j)-goiU,,)f = Opin-^^/^''+'^). 

i=l j=p+l 

(ii) With Assumptions 8 and 9 additionally, we have 
(3.3) V^0n-(3o)^N{O,^), 
where T, = K~^SK-^. 

4. Rank score test. The asymptotic normality of /3„ enables us to make 
large sample inference on the coefficients a^, bk and 70. However, the variance- 
covariance matrix S of Theorem 3.1 is not easy to estimate, as it involves 
the densities of e^j-. Likelihood-based tests such as those discussed in [10] 
are harder to develop for the quantile models than for conditional mean re- 
gression. The rank score test given by Gutenbrunner, Jureckova, Koenker 
and Portnoy [13] provided an attractive alternative to hypothesis testing 
on quantile regression coefficients by avoiding direct estimation of the error 
densities. In this section we extend the rank score test to the semiparamet- 
ric model (2.1). Fundamental theories of rank-based inference can be found 
in [12, 14, 28]. 
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As in [13], we assume Cjj to be i.i.d. in this section. Even so, the inclusion 
of a nonparametric component go in the model complicates the derivation 
of the limiting distribution of a rank score test statistic. Fortunately, the 
limiting distribution initially derived from linear quantile regression models 
remains valid when go is estimated by regression splines. 

4.1. Test statistic and asymptotic distribution. Suppose that Pqi is a q-di- 
mensional subset of /3o in (2.1), and is the N x q design matrix corre- 
sponding to Pqi. Let /3o2 and Xn2 be the remaining parameter vector and 
design matrix, respectively; then the global model (2.1) has a partitioned 
form 

(4.1) Yn = go{T) + Xj,f3o, + Xj^fB^, + e„. 

Suppose that we wish to test the null hypothesis Hq : f3Qi = 0, versus the al- 
ternative hypothesis Hi : (3qi ^ 0. Using the spline approximation described 
in Section 2, we further denote that 

Zn = {'^1,2: ^1,3^ ■ ■ ■ 5 ■^l^mi ) • • • ) ^n,2; ^n,3^ ■ ■ ■ ) ■^n,m„) ) 
Wn = iZj,X^,)^, 

= GnX^l, 

(/>o = (ao ,/3o2)'^, 

where oq is defined in Section 3 and Wn is the pseudodesign matrix when 
Hq is true. 

Let Vij be the column component of Vn that corresponds to the ith subject 
and the jth measurement. Similarly, we denote Wij as the column compo- 
nent of Wn that corresponds to the ith. subject and the jth measurement. 
The dimension of Wij is increasing in the dimension of the B-spline space. 

Now, let Sn = ri~^^'^Y.i=iY.]ll2V>T{Yij - wjj(i)n)vi,j, where = 

argmin,^ ^7=1 Ef=2 Pr{y^,i " <,'/') and K = n" V(l - t)X„i (/„ - Gn)Xji . 
We define the rank score test statistics as 

(4.2) Tn = SlV-^Sn. 

To obtain the asymptotic distribution of T^, we make additional assump- 
tions as follows. 

Assumption D1. The random errors e^j are independent of one an- 
other, and there exists a constant h such that hi^j = fe^ ^ (0) = b, for all i and 
j- 
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Assumption D2. ma-XijWwijW^ = Op{kn), E[maxij \\vij\\'^] < oo and 
Assumption D3. There exists a constant k such that sup^ j\ fy^ .^y. .^^^ . \ < 

K. 

In Assumption D3, we assume the conditional density of Yij to be uni- 
formly bounded. Assumption Dl holds if the e^j 's are i.i.d. Besides Assump- 
tions D1-D3, we have a slightly stronger assumption on the smoothness of 
the function go{t). We modify Assumption 1 in Section 3 as 

Assumption 1*. go(i) has bounded rth derivative for some r > 2. 

We now have the main theorem of this section, which allows, but is not 
limited to, the choice of kn ~ n^/('^^'+^) given in Theorem 3.1. Using <C bn 
to mean lim^^oo fln/^n = 0, we have the following result. 

Theorem 4.1. Under Assumptions 1*, 2-7 and D1-D3, if the number 
of knots kn for the B-spline space satisfies rL^I^"^ <^kn^ n}!^ with r > 2, we 
have Tn ^ xl as n — > oo. 

It is helpful to note that the limiting distribution of T„ is invariant over 
a wide range of choices for /c„ . Practically speaking, the result suggests that 
the number of knots is not a critical issue for the rank score tests. In our 
empirical investigations, the validity of the test appeared robust against 
violation of the i.i.d. assumption on Cij. 

4.2. An assessment with Monte Carlo simulations. We may assess the 
performance of the rank score test described in the preceding subsection by 
Monte Carlo. A set of models is selected with different error distributions, 
percentile levels r and model structures. In each model, the number of sub- 
jects is 200, the number of measurements from each subject is 10 and the 
percentile level is 0.5 or 0.95. The number of Monte Carlo samples used for 
estimating the Type I error is 10,000 in each scenario, so that the standard 
error of the estimate is around 0.2%. 

Model 1 A semiparametric longitudinal model with i.i.d. random errors. 
Vi,i = h{tij) + byij-i + Xi + ejj, where yij is the jth measurement from 
the ith. subject at time tij, h{t) =40t/(l -|-4t) is the intercept function, e^j 
are i.i.d. standard normal, Xi, independent of e^j, are normal with mean 
10 and standard deviation 1, and 6 is a parameter of interest. Here, h{t) 
has similar shape and scale as a reference centile curve of infant weight. 
Without loss of generality, we assume that, given i, the measurement times 
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Table 1 
Type I errors for Models 1-6 



# of knots 




T = 0.5 






T = 0.95 




Model 1 


Model 2 


Model 3 


Model 1 


Model 2 


Model 3 


1 


0.0524 


0.0535 


0.0531 


0.0508 


0.0572 


0.0523 


2 


0.0494 


0.0472 


0.0488 


0.0498 


0.0551 


0.0505 


3 


0.0513 


0.0485 


0.0479 


0.0495 


0.0550 


0.0500 




Model 4 


Model 5 


Model 6 


Model 4 


Model 5 


Model 6 


1 


0.0534 


0.0497 


0.0563 


0.0543 


0.0672 


0.0521 


2 


0.0506 


0.0462 


0.0505 


0.0504 


0.0662 


0.0517 


3 


0.0510 


0.0457 


0.0523 


0.0479 


0.0684 


0.0510 



ii,2) ■ • ■ ) ^i,io) ^-re order statistics of ten uniformly distributed random 
variables on the interval [0,1]. The intercept function h{t) is monotone in- 
creasing and is approximated by a cubic spline with internal knots (0.25, 0.5, 
0.75). Using the rank score test, we test the null hypothesis: Hq : 6 = 0, that 
is, the significance of the autoregressor is of interest. Models 2 and 3 have the 
same structure as Model 1, but assume that the distribution of ejj is t^/at^ 
and Xi/<7^2) respectively, where the variances of e^j's are standardized to 1. 

Models 4-6. A longitudinal model with non-i.i.d. errors. 

Vij = h{ti,j) + ^yij-i + + (1/2 + ?/ij_i/25)eij. 

For these models, the conditional quantile function is a linear function 
of Vi,j-i with the autoregressive coefficient as hr = b + Qr{eij)/25, where 
Qrisij) is the rth quantile of eij. With Models 4-6, we test the null hy- 
pothesis Hq -.br =0. 

In all the models, the intercept function h{t) is approximated by a re- 
gression cubic spline, which may vary with different knot selections. To see 
whether the rank score test is sensitive to the number of knots, we estimate 
h{t) in Models 1-6 with one or two or three uniformly spaced knots. Table 1 
gives the estimated Type I errors when r = 0.5 and r = 0.95 for the two sets 
of models. All of them are close to 0.05, the desired significance level. 

A look at the estimated power curves (not presented in this paper) con- 
firms that the performance of the test does not depend heavily on the choice 
of kn and is quite robust against modest deviations from the i.i.d. assump- 
tions on the errors. 

5. Example on infant weight. In this section we use the Finnish growth 
data [31] to demonstrate the use of the global model for subject screening 
and compare the results with the LMS-AR approach of Cole [6]. The data 
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we use here include measurements from a total of 1143 Finnish boys from 
birth to 2 years of age. 

5.1. Global model on infant weight. Monitoring the growth of weight 
during infancy is clinically important. For example, there is evidence that 
the rapid rates of weight gain during infancy could lead to obesity later in 
the childhood and may also be related to cardiovascular diseases later in 
life [37]. We use weight as the measurement for screening, but two prior 
measurements on weight and the current height are included in the model. 
Approximately, the measurements were taken monthly during this period of 
infancy, but few children were measured exactly on a monthly schedule. 

Let Wij and Hij be the j'th weight and height, respectively, of the ith 
subject at age tij. The global model (2.1) for the rth quantile can be further 
specified as 

2 

(5.1) Wij = gr{ti,j) + ^{ak,r + h,TDij^k)Wij_k + CrHij + eij. 

k=l 

Note that the error variables e^.j are r-dependent in the model. A set of 
quantiles with r € {0.03,0.1,0.25,0.5,0.75,0.9,0.97} is chosen for considera- 
tion. In the example, we use cubic splines with evenly spaced internal knots 
(0.5, 1, 1.5). 

Table 2 lists the estimates of ak^r, bk^r {k = 1,2) and Cr for the seven 
r's. The numbers in parentheses are the corresponding p-values based on 
the rank score test. They suggest that the second prior weight is not signif- 
icant at several r levels. Even when it is statistically significant, the small 
magnitudes of the coefficients mean that the contribution to the conditional 
quantiles would be small relative to that from the first prior weight. Height 
is a significant factor in the model. 

Also note that ai^r decreases, but bi^r and Cr increase with r. This means 
that the measurement time distance and height play larger roles for heavier 
boys than for lighter boys. 

A boy's growth status can be assessed by comparing his current weight 
with a set of conditional quantiles estimated from model (5.1). By boot- 
strapping subjects, we can also get confidence intervals on the conditional 
quantiles. 

The left panel in Figure 1 displays the weight growth path of one subject in 
the Finnish growth data. The growth path is plotted against the conventional 
(unconditional) centile curves of weight. We note that the weight of this 
subject had been growing along the 0.25th quantile but jumped to about 
the median level at age 0.61 year. A rapid weight gain could be a warning 
signal for obesity. Therefore, it is worthwhile to single him out early at age 
0.61. 
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Table 2 

Estimated parameters for infant weight 



T 





0.03 


0.10 


0.25 


0.50 


0.75 


0.90 


0.97 


ai,r 


0.787 


0.823 


0.810 


0.801 


0.737 


0.612 


0.482 




(0.000) 


(0.000) 


(0.000) 


(0.000) 


(0.000) 


(0.000) 


(0.000) 


6l,T 


0.163 


0.216 


0.279 


0.308 


0.332 


0.334 


0.346 




(0.000) 


(0.000) 


(0.000) 


(0.000) 


(0.000) 


(0.000) 


(0.000) 




0.012 


-0.015 


-0.016 


-0.032 


-0.048 


-0.056 


-0.037 




(0.496) 


(0.103) 


(0.015) 


(0.000) 


(0.000) 


(0.002) 


(0.295) 


62, T 


0.045 


0.020 


0.008 


0.023 


0.021 


0.043 


0.071 




(0.004) 


(0.083) 


(0.547) 


(0.080) 


(0.223) 


(0.012) 


(0.040) 


Ct 


0.038 


0.043 


0.051 


0.059 


0.080 


0.109 


0.138 




(0.000) 


(0.000) 


(0.000) 


(0.000) 


(0.000) 


(0.000) 


(0.000) 



The numbers in the parentheses are p-values. 
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Fig. 1. An example of conditional screening based on the global model. 



We apply model (5.1) to screen the weight at 0.61 year as the "current" 
measurement age. The right panel in Figure 1 shows the screening results 
from the global model, with the circle dots as the estimated conditional 
quantiles at age 0.61 from the global model. The boy's current weight is 
9.35 kg, which is higher than the 0.97th conditional quantile, suggesting 
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Fig. 2. A comparison between the global model and the LMS-AR model. 

that he might be overweight given his prior path and current height. The 
gray band represents the pointwise 90% confidence band for the conditional 
quantiles obtained from the bootstrap method. The subject's current height 
is higher than the upper bound of the 0.97th conditional quantile, which 
reinforced our conclusion that, from ages 0.46 to 0.61, the boy had been 
gaining weight at an excessive rate, even though his weight was still below 
the 0.75th quantile at age 0.61 based on the conventional growth chart. 

5.2. Comparison with the LMS-AR model. Cole [6] suggested using the 
LMS method to transform data to normality and then applying an autore- 
gressive model to the z-scores to estimate the conditional quantiles. To use 
the LMS-AR model, we interpolate the data to have measurements on fixed 
time intervals. More specifically, we first transform the infant weight using 
the LMS method with the effective dimensions 7, 10 and 7 for the functions 
L, M and S, respectively. We denote as zt the z-score at age t. An AR(2) 
model is then used at 0.61 year, with the first two lags at 0.46 year and 0.37 
year as autoregressors: 

(5.2) zo.ei = ao + aiZo.Ae + 02^0.37 + e, 

where e ~ A^(0, cr^). 

The estimated quantiles from the LMS-AR(2) model, together with those 
from the global model using the same covariates, are plotted in the left panel 
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of Figure 2. The two sets of conditional quantiles agree fairly well at the lower 
percentiles, although the former are slightly higher than the estimates from 
the global model. However, for the upper quantiles we observe considerably 
large differences between the two estimates: the estimates from the LMS- 
AR(2) model are much lower than the ones from the global model, and the 
difference widens as r increases. 

To understand the differences, we note that the LMS-AR approach im- 
plicitly assumes joint normality of the z-scores at three measurement times. 
The right panel of Figure 2 gives the normal QQ-plots of the z-scores for 
these four variables. While the LMS transformation does a pretty good job 
in achieving marginal normality at these three measurement times, it is clear 
that the distribution of the residuals deviates substantially from normality, 
and thus the joint normality is not achieved. The QQ-plot of the residu- 
als shows clearly a heavy upper tail, which explains why the LMS-AR(2) 
method underestimates the upper quantiles. 

6. Model diagnosis. To use the global model for subject screening, one 
has to have some faith in the adequacy of the global model. Does the model 
fit the data well? Can a more parsimonious model fit the data almost as 
well? These questions lead us to consider model assessment tools for the 
conditional model. 

Goodness-of-fit tests developed for linear quantile regression models by 
Koenker and Xiao [27] and He and Zhu [20] may be extended to semipara- 
metric models, but they are not made for model assessment over a set of 
quantiles. In this section, we propose an assessment tool tailored to our con- 
ditional models by trying to compare the empirical distribution of Y with 
the simulated distribution from the model. 

Because of the irregularity in measurement times, we would not have 
enough data for model assessment at any given time. However, we view Yij, 
the jth measurement of the ith subject, as a random variable of interest for 
any fixed j. We allow tij to be random, and consider the distribution of the 
jth measurement for any group of subjects. 

Given the quantile models, there is a nice way to generate Yij from its 
conditional distribution. Specifically, we draw U from the uniform distribu- 
tion on (0, 1), and then generate Y as the Uth quantile from the (estimated) 
model. A resulting random sample follows the model-based conditional dis- 
tribution. 

For a given j, we randomly draw a subject from those who have at least 
j measurements. Given this subject's y(ij_i) and X{tj-i), we draw Y{tj) 
from the quantile models as described above. When this process is repeated 
many times, we obtain a simulated sample for the jth measurement. The 
size of the simulated sample could be as large as we wish, but 5000 is used 
in our example below. If the model fits the data well, the distribution of the 
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simulated y{tj) should match the distribution of the observed jth measure- 
ments. It is a weak requirement that the two marginal distributions should 
match; instead we would expect a good match for any subgroup of subjects 
defined by the covariates in the model, which forms the basis of our model 
assessment tool. 

To have a quick idea of how well the empirical and the simulated distri- 
butions match, we may use the familiar QQ-plot. Another useful option is 
to look at \/n{T — t)/{t{1 — r))^/^ for a set of r's, where f is the propor- 
tion of the observed Yij that fall below the rth quantile of the simulated 
distribution, and n is the number of the observed Yij. 

For the infant weight example in the previous section, we display the QQ- 
plots in Figure 3 for j = 2, 4, 6 and 8. The reference line in each plot is y = x. 
These QQ-plots suggest a good match between the two distributions. The 
upper half of the figure plots the standardized f — r for 10 equally spaced 
r's from 0.05 to 0.95. No differences are worrisomely large, except that the 
match at j = 2 was not as satisfactory as the others. 

However, it is more meaningful to assess the agreement for subgroups of 
subjects. When all the subjects are included, we are only concerned with 
the marginal distributions, so even an unconditional model would do well. 

We now consider some subgroups based on birth weight or current height 
and compare three models: the unconditional model, the global model (5.1) 
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Fig. 3. Checking goodness-of-fit at the jth measurement, where N stands for the number 
of observations whose median measurement age (in year) is also displayed. 
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and a smaller global model obtained by dropping the second lagged term 
as well as height. Figure 4 gives the QQ-plots and the standardized (f — r) 
plots for two subgroups: Group 1 subjects whose birth weights are below 
the first quartile, and Group 2 subjects whose current heights are below the 
first quartile. Here, model adequacy is evidently different. The unconditional 
model, as expected, would overestimate the lower quantiles for Group 1 and 
Group 2. The smaller conditional model would fit the data better than the 
unconditional model, but the global model (5.1) does much better. The 
improvement due to lag-2 weight and due to height information is clearly 
demonstrated through the model assessment plots proposed in this section. 

By examining the model checking plots for various subgroups defined by 
covariates, we find that global model (5.1) is a good fit to the Finnish infant 
data for the purpose of conditional weight screening. Here, we choose not to 
specify a default for formations of subgroups; rather we leave it to the indi- 
vidual researcher to slice and detect. The model assessment ideas presented 
here are applicable to conditional quantile models in other applications, and 
formal lack-of-fit tests may be developed from these ideas. We defer them 
to future work. 
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7. Further discussion. It is clear that the conditional growth charts that 
we propose in this paper are not a set of charts in the usual sense. Screening 
based on the global model on longitudinal growth data needs to be done 
with customized charts for each subject. However, such customized charts 
can be generated easily given the estimated models, and all one needs is 
to input the information on a given subject under screening, making the 
clinical use of the conditional growth charts really practical on any desktop. 

In this section we discuss some limitations of and generalizations to the 
model (2.1). We hope that our work will stimulate future research on longi- 
tudinal growth charts. 

7.1. Generalization. Fixed measurement intervals are usually assumed 
in an autoregressive model. By assuming that the AR coefficients are linear 
in the time spacings, the global model (2.1) has some flexibility to deal with 
irregular measurement times. We can further generalize the model as 

p 

(7.1) Yij = goiUj) + hk{Dij^k)Yij-k + Xj^jlo + eij, 

k=l 

where the AR coefficients are some unknown functions of measurement time 
spacings. In fact, for small measurement time spacings, + bkDij^k works 
reasonably well as a linear approximation to the general function hk{Di j^k). 
The estimate of hk{Di^j^k) in model (7.1) can be obtained by either a para- 
metric or a nonparametric approach. For a parametric approach, we may 
specify a parametric coefficient function, such as + h^Di j^k + ^kDf ^ ^. To 
choose between + bkDij^k + ^kDf j ^. and -I- bi^Dij^k, we can simply test 
the hypothesis Hq : = 0. For a nonparametric approach, one may use the 
varying-coefficient models explored by Kim [25] or the dynamic models of 
Cai and Xu [2]. 

Another restrictive feature of the global model is that the AR coefficients 
do not depend on the measurement times, but only on their spacings. An 
alternative model can be set up as 

p 

(7.2) Yij = go{tij) + Yi^-kitiJ-k) + bk{tij-k)Dij^k]Yij-k + X^jlo + Cij, 

k=l 

where we allow the parameters and bk to be functions of prior mea- 
surement times tij-k- Time series models studied by Xu [44] may also be 
considered. In practice, the flexibility of the models becomes an advantage 
only when the sample size is large. 

7.2. Localization. If all the subjects have measurements on a series of 
fixed time points, say {to,ti,t2, ■ ■ ■ ,tp), then and bk in the global model 
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are no longer identifiable. To model the conditional quantiles of the mea- 
surement at time to given the prior path and other covariates, one may 
use 

p 

(7.3) Yt^^i = (3o + J2PkYt^,r + l'^X, + e^, i = l,...,n, 

k=l 

where Yt^i is the measurement of the ith subject at time t, (ti,t2, ■ ■ ■ ,tp) are p 
prespecified earlier time points, which are not necessarily evenly spaced, Xi is 
an /-dimensional covariate and Ci are i.i.d. random variables with zero rth 
quantile. It follows that the conditional rth quantile of Yt^^i is a linear com- 
bination of p prior measurements ^2,*' • • • ' ^p,«> the other covariate 
Xi. 

7.3. Conditional versus marginal growth charts. We have focused on the 
conditional growth charts in the present paper, because conditional charts 
are highly valuable in evaluating recent growth status of subjects. However, 
we do not suggest that we should do away with marginal growth charts, 
which are most commonly used by health professionals today. 

The conditional growth charts may be unsuccessful in screening out sub- 
jects with gradual but persistent slowdown in growth. Further studies are 
needed to address issues related to the "longitudinal" use of growth charts 
over time. If a subject consistently falls below the 25th conditional per- 
centile for several time periods, it might signal a growth problem. At the 
present time, the marginal growth charts may be used in conjunction with 
conditional charts. 

7.4. Monotonicity. The growth charts from model (2.1) require estima- 
tion of a series of conditional quantile functions. It is not guaranteed that 
the estimated quantiles are monotone in r. For a discussion of constrained 
quantile regression that avoids this crossing problem, see [16]. 

8. Technical proofs. 

8.1. Proof of Theorem 3.1. First we write model (2.1) in matrix form as 
(8.1) Yn = go{T) + Xjf3o + en, 

where T = {tij)Nxl = (^1,2, . . . ,tl,mi,t2,2, ■ ■ ■ ,t2,m2i ■ ■ ■ 5*n,2, . . . ,*n,m„)''', = 

{Yij)Ny.i, Xn = {Xlj)Nx{i+2), dn = (eijOiVxi and (3q = (00,60,70" Y . The in- 
dices (i,j) in the vectors Yn and e„ and the matrix Xn are arranged in the 
same way as in T. Here N = Y^^=i{''^i ~ 1) is the total number of observa- 
tions used for the estimation. Note that rui are uniformly bounded for all 
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n, 0{N) = 0{n). Recall that Z„ = (7ri^2, 7ri,3, • • • , tTij, • • • , vr,,, rn.. iTr^i, and 



Nxkn 



Gn = Zn{Zj Zn)Zj . We further define 

Note that z„ 7^ Z„ in our notation; Zn is the matrix of the spline basis func- 
tions, and Zn is the whole design matrix after normalization and orthog- 
onalization. We split z„ into two parts: the matrix Zni is the normalized 
design matrix of the spline bases, and Zn2 is the normalized design matrix of 
the linear components that is orthogonal to Zni- Also, we denote Zni,ij and 
Zn2,ii as the column components associated with the ith subject and the jth 
measurement. Similar partitions and standardizations were used in [21]. 
With this notation, model (8.1) can be further written as 



Yn = Znl6i{0LQ,l3Q) + Zn262{OLO-,(io) + Rn + en 
= ZnO{cXo, (3q) + Rn + en- 



The orthogonality between Zni and Zn2 will not affect the estimation of 
Qq and Pq, but could simplify our proof. Using this orthogonalized model 
format, let 

C'ij(^) — Pri^ij Zn^ijO Rnij) Pri^iJ Rnij) ~\~ ^n^ijOlpl^eij^ , 
Aj (^1, ^2) = Pri^iJ — Znl,ijOi — Zn2,ij02 — Rnij) 

— Prie-ij — Znl,ijOi + Rnij) + ^n.2,ij^2'0(ei,j ) ■ 

The proof of Theorem 3.1 relies on the following three lemmas. 

Lemma 8.1. Under Assumptions 1-7, for any sequence {Ln} such that 
1 < Ln < for some 60 G (0, r/10), we have 



i.3) 



sup /c„ ^ 

9\\<Lr,kl/^ 



n nii 



i=lj=2 



Op(l), 



and for any Mi > and M2 > 0, we have 



(8.4) sup 

||6»i||<Afifcy^;||6»2||<Af2 



5:^(A,,(^i,e2)-MA,,(^i,e2)]) 

i=lj=2 



:Op(l). 



Proof. The uniform convergence of (8.3) and (8.4) can be proven fol- 
lowing arguments similar to those used in [19]. We skip the details to save 
space. □ 
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Lemma 8.2. Under the same conditions of Lemma 8.1, we have, as 
n — > oo, 



\\e\=i " 



Y^Y.^E[C,ALnkl/^0)] 

i=lj=2 

— L^kJ Zn^ijOip{eij)) 



>n-i. 



Proof. First, under the condition of \\9\\ = 1 and d„ < W5, it follows 



from Assumption 5 that E[Cij{9)] = E[J_j^ 

where Bn is a bounded sequence. Let dn = maxjj(L„,A;' 
we have 



{bijs + Ens"^) ds], 



1/2|| I 
n \\Zn,ij\ 



+ \R 



mj I ; ) 



(8.6) 



i=lj=2 



i=lj=2 



hi 
2 



E 



i=l j=2 



i=lj=2 



n mi 



+ Lnk^ ^^"^ ^ ^ E{bijRnijZn,ijO) 
i=lj=2 



+ 0(1). 



1/2 



To see the last step of (8.6), we first note that = 0{n ^^'^k}/'^), 

\\e\\ = 1, and thus kn^^'^Y,1^^J2Y=2E\^n,ijd\ < oo- Meanwhile, d^L^ < d„ 
Op(l). It follows from the dominated convergence theorem that E[d'^Lnkn 
J27=i YL"j=2 \^n,vfi\\ = 0(1), and therefore the last step of (8.6) holds. 

Since hi^j is uniformly bounded away from zero and infinity due to As- 
sumption 5, it is easy to see that Y^=\ bounded away 
from and from infinity, uniformly in ||0|| = 1. On the other hand, we have 



-1/2 ■ 



Y^=\Y^=2 ^{bi,jRnijZn,ij(^) = 0{Ln)- Therefore, we conclude that 



k-^Y.T.E[C,,{Lnkl/H)]>cLl 

i=lj=2 
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for some c > 0, which imphes that we only need to show 

n mi 

(8.8) LnK^I^ E ^n,^je^l^{e^,J) = OjXLl) 

i=lj=2 

to prove Lemma 8.2. 

By Theorem 4 of [3], the eigenvalues of knZ^ Zn/n are bounded away from 
zero and infinity when n is sufficiently large. Assumption 4 implies that the 
eigenvalues of n~^ipT-{en)ipT{en)'^ are also bounded. Based on these facts, if 
we denote Vr(en) = {^pT{el^l), . . . ,ipr{en,mjV , then Vr (en)^V'r(e„) = 
diag{Tpr{ei)'^ 'ipriei)) ■ Let Aj be the largest eigenvalue of V'T(ei)''''(/'r(ei); supj Aj 
Op{l) by Assumption 4. It thus follows that 

n mi 

sup LnK^^'^YY Zn,ij9lpT{eij) 
m\=l i=lj=2 



(8.9) 



sup L„/c„^/2[6''^zJ'^^(en)^rC 
l|e||=i 



< sup supA,L„A:-i/2[0T^T^^^]i/2 
||e||=i i 



^supA,jL„A;-i/2 = Op(L„fe-^/2)_ 
Therefore, (8.8) holds, and the proof of Lemma 8.2 is complete. □ 



Lemma 8.3. Under Assumptions 1-9, we have 

n mi 
i=lj=2 



(8.10) sup 

\\9i\\<Lk'J^;\m<M 



Op{l). 



Proof. First of all, we note that 

n m.i 

^^i?[Aj(ei,e2)] - \n-^ejKne2 

i=lj=2 



sup 

9i||<Lfcy2; 
I16'2||<M 



sup 

h\\<Lkl/^; 



i=lj=2 ^ 



-Z„l^ij0l—Z„2,ijd2—Rn 
^7il ,i 7 ^1 Rnii 



i>{eij + s)ds 



i\\<M 



(8.11) 



sup \^[E{e2 Z„2BnZn202) " 6*2 2n2^n2n26'2] + ^(6*2 Zn2rr. 
\92\\<M 

+ 0(1), 
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where r„ is the n x 1 vector with elements (bijRnij)- 
Assumption 8 imphes that 

(8.12) sup \E{9j Z^^BnZn202) " ^2^42^n^n2^2 1 = Op(l). 

\\S2\\<M 

On the other hand, sup||g2|j<A/ ^2~^^2'"n ^^P\\02\\<M ^2 ^n2'^ni and it fol- 

lows that 

(8.13) sup E[eJz^^rn] = o{l). 

||6»2||<M 

Combining (8.11)-(8.13), we arrive at (8.10). Lemma 8.3 is thus proven. □ 

With Lemmas 8.1 to 8.3, the proof of Theorem 3.1 can be outlined as 
follows: 

Proof of Theorem 3.1. Model (8.1) can be reconstructed as 

Yn = Znl0i{aQ,(3Q) + Zn202{oL0,PQ) + Rn + 
= Zn6{ao, (3q) +Rn + en- 

We define 

= (^nli ^n2)^ 

(8-14) ^ „ ^ 

= (6'i(q!„,/3„) - di{ao,/3Q),e2{an,i3n) - ^2(0:0, /3o)) • 

The objective function (2.2) can be written as 

n rrii 

(8.15) ^ ] y ] Pri^iJ ~ Zn,ij0n ~ Rnij)^ 
i=lj=2 

and we have the fact that /3„ minimizes (2.2) and On minimizes (8.15). 
Directly from the definition of On, we have 

(8.16) 0n - M = \K^^^0n2\\ = Op(n"^/2||0~„,2||), 

— 1/2 

which implies that \\0n2\\ = Op{kn ) is a sufficient condition for (3.1). On 
the other hand, 

n rui 



1=1 ]=2 



2 n mi 



<-EEK.("«-«o))' + 2t^|fe 

^.17) 



-2r 



1 ■ o 

1=1 3=2 



< 2[n-'\\0ni\\' + n"'\\H-'k'J'zJXn0n - /3o)f ] + 2W^l^n"' 
= Opin-'\\0nif) + 0(0- Poll) +Oikn^''). 
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The last step of (8.17) is due to the fact that n '^\\GnXn\\'^ < n ^\\Xn\\'^ < oo 
under Assumption 7. Combining (8.16) and (8.17), we conclude that it is 

~ —1/2 

sufficient to show ||0n|| = Op(A;„ ). 

According to Lemmas 8.1 and 8.2, for any e there exists such that 

{n TTii n rrii "\ 

1/2 X! P-ri^iJ ~ ^'n,ijd - Rnij) >YY1 P-r^^hj ~ ^riij) > >l-e. 
\\e\\>Leki'^ i=lj=2 j=lj=2 J 

Since 6n minimizes J27=iJ2Y=2 Pr{^i,j ~ Zn,ijS — Rnij) over the space R^", 

we have P(||0„|| < L^kn"^) > 1 — e, and thus \\6n\\ = Op{kn'^). 

With the consistency of the global model, we can take a further step to 
show the asymptotic normality of the coefficient estimate 

We denote by z„2 the submatrix of Zn2 corresponding to the ith subject; in 

other words, Zn2 = {z^nF ^ ^2^' • • • ' 4?'^)^- Let O^^ = nK'^ Er=i ^^^(ei)- 
Due to Assumptions 8 and 9, 6^2 asymptotically normally distributed with 
asymptotic variance-covariance matrix K~^SK~^ . Since 6n2 = n^^'^Pn^ if we 
can show that 



(8.18) \\9n2-0*n2\\=Op{l), 

then (3.3) holds. 

Due to the definition of and the consistency of On, we know that 
P{ei2 <M)^1 and P(||^„,i|| < Lki^'^) 1 for any L > and M > 0. Let 

Aj(^2,^2) = Pri^iJ — Znl,ijOnl — Zn2,ij02 — Rnij) 

Pri^^iJ Znl^ijOnl Zn2,ij02 Rnij)- 

It follows from Lemma 8.4 that, for any given (5 > 0, 

n nii 

EE{A,,(^2,C2) - Zn2,ij^{e,,,){e2 " ^2) 



sup 

\\e2-e*^J<i 



i=l j=2 



^[Aj(^2,C2)]} 



Op{l). 



Furthermore, by Lemma 8.3 we have 

n rrii 

Y.Y.yD^,{02,el2)] 

i=lj=2 

n 



sup 

\e2-e*^J<& 



+ (^2-^:2)^E4TV'(e^ 



n 



-1 1/3T 
2^2 



Kn02 + n ^l02n^n'yn2 



i=l 



(8.19) = sup 

\S2~S*J<5 



Y.Y.iD.A02,e:2)] 

i=lj=2 
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sup 



1^2- 
Op{l). 



\\<& 



+ n ^[{62 - 0^2)^ KnOl - \6jKnO2 + \G2J1KnOn2] 



i=lj=2 



For sufficiently large n, n ^{62- 6^2)^ Kn{02 - 6*^2) > when \\62 - 6**2!! > 
5. Then (8.19) implies that 



lim P\ 



inf 



^.20) 



n rui 
^2\\>\=ij=2 



n mi 



60 — R 



^ ^ Prifii^i ^nlyij^nl ^nl,ij02 Rnij] 
i=lj=2 



1. 



Since 9ni minimizes X^Li Ej^2 Prieij - Zni,ij6i - Zni,ij92 - Rnij) over R'+2, 

(8.20) implies that for any 6, P{\\en2 - 6'n2ll > "5) ^ 0, that is, 11^2 - 6*^211 = 
Op(l). The proof of Theorem 3.1 is hence complete. □ 

8.2. Proof of Theorem 4.1. Recall that 5„ = n'^l'^ ELi E?^2 V'r (>ij - 



Since the t'ij 's are the least squares residuals from regressing Xni 
on Wn, they differ from vf^^ = Xni,ij - E{Xni,ij\'Wn,ij] by Op{k:^''). It is easy 
to show that the limiting distribution of T„ = S^V~^Sn will not change if 

Vi^j are replaced by vf^; the latter enjoy intersubject independence and are 
often easier to handle mathematically. To simplify notation, we will simply 
prove the results in this section by assuming that the fjj's are intersubject 
independent. 

First, we give the following two lemmas. 



Lemma 8.4. Let Uij{(j), 0o) = 'PriYij - Wij(p)vij - (priYij - wljMvij - 
Eipr {Yij — wjjcl))vij + Eip-r {Yij — wJjcl)o)vij . Under the assumptions of The- 
orem 4.1, for any L > 0, we have 

n rrii 



i.21) 



sup 

-0o||<L(fe„/n)V2 



i=l j=2 



OJn^l'kT{Hn)f") 



An-'"). 



Proof. Lemma 8.4 can be viewed as a special case of He and Shao [17], 
which considered the asymptotic behavior of M-estimators with increasing 
dimension. 
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Since (frit) = 1 — Tl^t<o} is a constant function with a jump at i = 0, we 
have 



^.22) 



<M\vLiri 



{ I y,,, -wjM < . {4>-<Po ) 1} ■ 



Under Assumption D3, 



^.23) EJ2 



i=l 



i=2 



, \ n mi 

^ * ^ i=l,=2 



which, together with Assumption D2, imphes 
(8.24) A 



sup E4 

U-M\<L{k„/ny/^ i=l 



'Yuij{(j),(j)o) 

i=2 



0((?iA:„)i/2). 



Moreover, since v'T(i) is bounded and maxjj = it follows from 

Lemma 2.2 of [17] that 



^.25) Bn= sup 



J =2 



Finally, combining (8.23) and Assumption D2, we note that condition (CI) 
of Lemma 3.3 of [17] is satisfied, and (8.21) holds consequently, where the 
right-hand side comes from Op((ln(n)A;„)^/^(n~^ + aI/"^ + bU'^)). □ 

Lemma 8.5. Under the assumptions of Theorem 4.1, 

n mi 



^.26) 



sup n 

-9io||<i(fe„/n)l/2 



-1/2 



j=lj=2 



Eipr{Yij - wJj(j)o)ViJ 



oil). 



Proof. Under the condition ||(/' — (^oll < Likn/nY/"^ , we first expand the 
conditional mean E[(pr{Yi^j — wJj(l))vi,j\vi,j,Wij] around (po for each i and j. 
Note that, for a real-valued random variable u E[ipT-{u)] = r — -Fu(O), where 
Fu is the distribution of u, we have 

E{[ipr{Yij - wjj(j)) - (PriYij - wlj4>o)]Vij\ViJ,Wij} 

(8.27) =fe,JO)v,,jwl^{(^-(^o) 
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The remainder term in (8.27) is o{n ^) because — 4>q\\ < L{kn/n)^^'^ and 
kn = o(n^/^). It follows from (8.27) that the left-hand side of (8.26) is 



sup n 

U-<f>o\\<L{kn/ny/^ 



-1/2 



n nii 



i=lj=2 



i.28) = sup n~^/^\\bEVjWn{ 

ll<^-<^oil<L(fc„/n)i/2 



+ 0(1), 



where b = fei^{0) and c = /e^^, (0). Note that Vn and Wn, are orthogonal to 
each other, that is, VjWn = 0, the first term in (8.28) is zero. Moreover, 

sup \E{(I) - ct)oyWn[diag{vi,j)]Wj{^ - M\ 



<L^{kn/n)E 



rnaxf ||yV„,| 



0{kl). 



The last equation holds due to Assumption D2 and the fact that £'||Wn|P 
0{nkn). When k^ = o{n^/^), 

(8.29) n-^l^cE{(t>-(t>G)^Wnd\B.g{vij)Wl{<t)-(t>G) = o{\). 

Thus, (8.26) follows from (8.28) and (8.29). □ 



Proof of Theorem 4.1. In Theorem 3.1 we have shown the consis- 
tency of given kn ~ v}^^'^'''~^'^\ which, however, is not the necessary condi- 
tion for the consistency. In fact, for any fc„ = o{'n}^^), the consistency of (pn 
holds at the convergence rate {kn/nY/"^ , that is, ||0„ — 4>q\\ = Op{n~^/'^kn'^). 
Therefore, we can derive from Lemmas 8.4 and 8.5 that, when kn = o(n^/^). 



n 



-1/2 



^.30) 



i=l j=2 



- iPr{Yij - wJj(l)o)Vij} 



Opil). 



Denote 5* = J2i=i Ej^2[v'r(eij>ij]; then the summands J2f=2[^Tieij)vij] 
are independent of each other and have mean zero. Due to the between- 
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subject independence, 

n / mi \ n 

i=l \j=2 ) i=\ 

n 

= r(l - r) ^ V7V. = t(1 - r)X„i(/„ - 

where Vj is the rrij x rrii submatrix of Vn corresponding to the zth subject, 
and (friei) is the same as defined at the beginning of Section 3. In the i.i.d. 
case, Var(99^(ej)) = r(l - t)/^^. Let K = - r)X„i(/„ - it 

follows from the CLT that 

n mi 

(8.31) n-'^^Y.T.i^r{ei,j)vij] isA7V(0,K). 

i=lj=2 

Combining (8.30) and (8.31), it is clear that all we need to show Theo- 
rem 4.1 is 



^.32) n 



-1/2 



i=lj=2 



Op{l). 



However, it is hard to show (8.32) directly for some kn, more specifically, for 
^i/4r <^ /^^ ^ n^/"^^ including /c„ ~ n^^'^^^^. Instead, we take an intermediate 
step. 

Let Kn be the set of knots used in our global model estimation. We assume 
that the dimension of Kn ~ the same proof goes through for other 

values of Kn- Under Hq, let Yij — 'wJj4>o = Rnij + eij, where Rnij is the bias 
from the B-spline approximation using the set of knots Kn- The bias will 
go to zero with supj ^ \ Rnij\ < W^k'^ , where W3 is defined in Assumption 1. 
By adding more knots into Kn, we have a new set of knots Kn with its 
dimension, denoted as kn, approximately n", where l/(2r) < a < 1/4. We 
also assume that the knots are added in such a way that the extended set of 
knots is quasiuniform. Using the new set of knots, we define Wn, Wij, Vn, 
Vij, (j), (t>o and Rnij the same way as >V„, Wij, V„, Vij, 0„, (po and Rnij in 
the original setting. 

Note that Vij is the residual of a spline estimate for the (i,j)th ele- 
ment of £^[X„i|Wn], using Kn as the set of knots; while Vij is the residual 
of a spline estimate for the (z,j)th element of £'[X„i|VV„], using Kn as 
the set of knots. Therefore, by the same arguments used for Theorem 3.1, 
^J2iJ2j = Op{kn'^^) + Op{kn/n). In the meantime, by definition, 

Vn is orthogonal to both Wn and Wn- Based on the facts above, following the 
similar arguments used for (8.30), we can show that, for n^^'^^ <^kn^ n^^^ 
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and r > 2, 
(8.33) 

and 

(8.34) n"^/2 



n rrii 



i=lj=2 



n nii 



i=lj=2 



Op{l) 



Op{l). 



Moreover, noting that, given Vij, ELpriYij — wJj(j)Q) = Eipr{Rnij + ejj) = 

Op{l). 



0{Rnij), we also have 
(8.35) 



i=lj=2 



Combining (8.33)-(8.35), we have 

n mi 



i.36) n-V2 



i=l j=2 



Under Assumption D3, and the fact that |(^T-(eij + -Rmj) — V'r(ejj)| < /||g. .|<;|^ 
we have 



n 



^.37) 



i=l j=2 

n rrii 

< 2qn-^/^EYT. \^nij\\Kj\\ = Oik-'n^/^) = o(l). 

i=lj=2 



The result (8.32) follows immediately from (8.30) and (8.37). The proof of 
Theorem 4.1 is thus complete. □ 
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